%%%%%%%%%%%%%%%%%%
% This is a TEMPLATE script for running the manual intervention OSL recommended preprocessing
% pipeline on Elekta-Neuromag data (a very similar pipeline will work on
% CTF data as well). 
%
% You'll need to do alter (at the very least) the settings in:
% datadir, workingdir, fif_files, spm_files, structurals_files

%%%%%%%%%%%%%%%%%%
%% SETUP THE MATLAB PATHS
% make sure that fieldtrip and spm are not in your matlab path

global OSLDIR;
    
tilde='/home/mwoolrich';
%tilde='/Users/woolrich';
osldir=[tilde '/homedir/matlab/osl1.3.1'];    

addpath(osldir);
osl_startup(osldir);

%%%%%%%%%%%%%%%%%%
%% SPECIFY DIRS FOR THIS ANALYSIS

% directory where the data is:
datadir=[tilde '/homedir/vols_data/murphy_project/data'];

% this is the directory the analysis files will be stored in:
workingdir=[datadir '/spm_files']; 
cmd = ['mkdir ' workingdir]; unix(cmd); % make dir to put the results in
 
%%%%%%%%%%%%%%%%%%
%% Set up the list of subjects and their structural scans for the analysis 

clear fif_files spm_files_basenames structural_files

% Specify a list of the existing raw fif files for subjects
% Note that here we only have 1 subject, but more generally there would be
% more than one, e.g.:
% fif_files{1}=[testdir '/fifs/sub1_face_sss']; 
% fif_files{2}=[testdir '/fifs/sub2_face_sss']; 
% etc...
raw_files{1}=[datadir '/fifs/hntface']; 

% Specify a list of the base fif file for use post Maxfilter
% Note that here we only have 1 subject, but more generally there would be
% more than one, e.g.:
% fif_files{1}=[testdir '/fifs/sub1_face_sss']; 
% fif_files{2}=[testdir '/fifs/sub2_face_sss']; 
% etc...
fif_files{1}=[datadir '/fifs/hntface']; 

% Setup a list of SPM MEEG object raw base file names to be created (in workingdir), in the same order as spm_files and fif_files:
% Note that here we only have 1 subject, but more generally there would be
% more than one
spm_files_basenames{1}=['spm8_meg1_new.mat'];

% Setup a list of existing structural files, in the same order as spm_files and fif_files:
% Note that here we only have 1 subject, but more generally there would be
% more than one.
structural_files{1}=['']; % leave empty if no .nii structural file available

%%%%%%%%%%%%%%%%%%
%% Create 1st round of unMaxfiltered data
% NOTE: Maxfilter does not work if the output fif file already exists!!!
for subnum=1:length(fif_files)
    disp(['%%%%%%%%%%%%%%%%%%%%%%%  SUBNUM = ' num2str(subnum) '  %%%%%%%%%%%%%%%%%%%%%%%'])
    Smf=[];
    Smf.fif=raw_files{subnum};
    Smf.fif_out=[fif_files{subnum} '_nosss'];
    Smf.logfile=1;
    Smf.downsample_factor=4;
    Smf.nosss=1;
    fif_sss=osl_call_maxfilter(Smf);
end

%%%%%%%%%%%%%%%%%%%
%% Convert to SPM

S2=[];
for subnum=1:length(fif_files), % iterates over subjects
    disp(['%%%%%%%%%%%%%%%%%%%%%%%  SUBNUM = ' num2str(subnum) '  %%%%%%%%%%%%%%%%%%%%%%%'])
    S2.spm_file=[workingdir '/spm8_' fif_files{subnum} '_nosss'];
    S2.fif_file=[fif_files{subnum} '_nosss.fif'];
    [D spm_files{subnum}] = osl_convert_script(S2);
end
close all

%%%%%%%%%%%%%%%%%%%
%% Look for scanner artefacts (bad channels) using OSLVIEW

for subnum=1:length(fif_files)
    D=spm_eeg_load(spm_files{subnum});
    oslview(D);
    waitfor(gca)
end

%%%%%%%%%%%%%%%%%%%
%% Re-Maxfilter data using bad channel information

for subnum=1:length(fif_files)
    disp(['%%%%%%%%%%%%%%%%%%%%%%%  SUBNUM = ' num2str(subnum) '  %%%%%%%%%%%%%%%%%%%%%%%'])
    Smf=[];
    Smf.fif=raw_files{subnum};
    Smf.fif_out=[fif_files{subnum} '_sss1'];                %set filename depending on which maxfilter option used
    Smf.logfile=1;
    if(~Smf.movement_compensation)
        Smf.downsample_factor=4;                             %if running movement compensation downsample has to be swtiched off
    end;
    Smf.nosss=0;
    Smf.spmfile=[workingdir '/' spm_files{subnum}];
    Smf.movement_compensation=1; %can only run with sss AND downsampling has to be switched off
    fif_sss=osl_call_maxfilter(Smf);
end

%%%%%%%%%%%%%%%%%%%
%% Convert to SPM

% set filenames used in following step
for subnum = 1:length(fif_files), % iterates over subjects
    spm_files{subnum}=[workingdir '/' spm_files_basenames{subnum}];
end

S2=[];
for subnum=1:length(fif_files), % iterates over subjects
    disp(['%%%%%%%%%%%%%%%%%%%%%%%  SUBNUM = ' num2str(subnum) '  %%%%%%%%%%%%%%%%%%%%%%%'])
    S2.spm_file=spm_files{subnum};
    S2.fif_file=[fif_files{subnum} '_sss1.fif'];
    [D spm_files{subnum}] = osl_convert_script(S2);
end
close all

%%%%%%%%%%%%%%%%%%%
%% Downsample with SPM (particularly important if movement compensation 
% used - but worth doing anyway)

%set filenames used in following steps
for subnum = 1:length(fif_files), % iterates over subjects
    spm_files{subnum}=[workingdir '/' spm_files_basenames{subnum}];
    %spm_files{subnum}=[workingdir '/d' spm_files_basenames{subnum}];
end

S=[];
for subnum=1:length(fif_files), % iterates over subjects
    disp(['%%%%%%%%%%%%%%%%%%%%%%%  SUBNUM = ' num2str(subnum) '  %%%%%%%%%%%%%%%%%%%%%%%'])
    S.D=spm_files{subnum};
    S.fsample_new = 250;
    D = spm_eeg_downsample (S);    
end
close all

%%%%%%%%%%%%%%%%%%%
%% Visual Inspection with OSLView (to find really bad artefacts in continuous data BEFORE Artefact removal with Africa)

%set filenames used in following step
for subnum = 1:length(fif_files), % iterates over subjects
    spm_files{subnum}=[workingdir '/' spm_files_basenames{subnum}];
    %spm_files{subnum}=[workingdir '/d' spm_files_basenames{subnum}];
end


for subnum=1:length(fif_files)
    D=spm_eeg_load(spm_files{subnum});
    oslview(D);
    waitfor(gca)
end

% At this point in the process, if desired you could perform semi-automated 
% ICA denoising using AFRICA. 
if(1) 
    %%%%%%%%%%%%%%%%%%%
    %% Perform AfRICA

    %set filenames used in following step
    for subnum = 1:length(fif_files), % iterates over subjects
        spm_files{subnum}=[workingdir '/fd' spm_files_basenames{subnum}];
        %spm_files{subnum}=[workingdir '/d' spm_files_basenames{subnum}];
    end

    for subnum= 1:length(spm_files),
        disp(['%%%%%%%%%%%%%%%%%%%%%%%  SUBNUM = ' num2str(subnum) '  %%%%%%%%%%%%%%%%%%%%%%%'])

        S=[];
        
        % set identification settings
        S.ident.eog_chan=312;                      % set to the right channel (EOG or eyetracker channel) 
        S.ident.do_blinks= 1;                      % does automatic blink detection           
        S.ident.do_ecg= 1;                         % does automatic ECG detection  
        S.ident.do_mains= 1;                       % does automatic mains artefact detection  
        S.ident.do_kurt=-1; 
        S.ident.ecg_corr_thresh=0.15;
        S.ident.manual_approval=1;
        S.ident.func=@identify_artefactual_components_manual;
        S.used_maxfilter = 1;                % indicates if maxfilter has been used            
        
        spm_file=[opt.dirname '/' spm_files_basenames{subnum}];
        S.fname=spm_files{subnum};                    
        S.logfile = 1;                                 
        S.ica_file = [S.fname '_preproc_ica_results'];           
        S.do_plots=1;                                          
        S.just_ica=0;       % set to 1 if want to run ICA only and not other two stages as well  

        S.to_do=[1 1 1];    %  here define which ICA stage to be run 1 = ICA only, 2 = component classification, 3 = component removal 
                 
        [spm_files_new{subnum}]=osl_africa(S);

        close all    
    end
end;

%%%%%%%%%%%%%%%%%%%
%% High Pass Filter
% To get rid of low frequency fluctuations in the data

%set filenames used in following step
for subnum = 1:length(fif_files), % iterates over subjects
    spm_files{subnum}=[workingdir '/A' spm_files_basenames{subnum}];
    %spm_files{subnum}=[workingdir '/Ad' spm_files_basenames{subnum}];
end

for subnum=1:length(spm_files),
    disp(['%%%%%%%%%%%%%%%%%%%%%%%  SUBNUM = ' num2str(subnum) '  %%%%%%%%%%%%%%%%%%%%%%%'])
    S=[];
    S.D = spm_files{subnum};
    S.filter.band='high';
    S.filter.PHz=1;
    S.filter.dir='twopass';
    D = spm_eeg_filter_v2(S);
end

%%%%%%%%%%%%%%%%%%%
%% Visual Inspection with OSLView (again to make final data clean before it goes into ICA, get rid of blinks, muscle artefacts etc that AFRICA didn't take care of)  

%set filenames used in following step
for subnum = 1:length(fif_files), % iterates over subjects
    spm_files{subnum}=[workingdir '/fA' spm_files_basenames{subnum}];
    %spm_files{subnum}=[workingdir '/fAd' spm_files_basenames{subnum}];
end

for subnum=1:length(fif_files)
    D=spm_eeg_load(spm_files{subnum});
    oslview(D);
    waitfor(gca)
end

%%%%%%%%%%%%%%%%%%%
%% OPTIONAL: DO EPOCHING (if epoch-based task data)
%% DO preliminary epoching for the purpose of finding outliers
% This is not the final epoching. Instead this sets up the epoch
% definitions, and performs a temporary epoching for the purpose of doing
% semi-automated outlier trial rejection (before running the fully 
% automated OAT).
%
% The epoch definitions and the continuous data will be kept and
% passed into OAT. This is so that things like temporal filtering (which is
% dones as part of OAT) can be done on the continuous data, before the data
% is epoched inside OAT.
%
% Note that this will also remove those trials that overlap with the bad
% epochs identified using oslview. 
%
% Here the epochs are set to be from -400ms to +500ms relative to triggers
% in the MEG data, We also specify the trigger values for each of the 4
% epoch types of interest (motorcycle images, neutral faces, fearful faces,
% happy faces). 

%set filenames used in following step
for subnum = 1:length(fif_files), % iterates over subjects
    spm_files{subnum}=[workingdir '/fA' spm_files_basenames{subnum}];
    %spm_files{subnum}=[workingdir '/fAd' spm_files_basenames{subnum}];
end

for i=1:length(spm_files), % iterates over subjects

    %%%%
    % define the trials we want from the event information
    S2 = [];
    S2.D = spm_files{i};
    D_continuous=spm_eeg_load(S2.D);
    
    S2.pretrig = -2000; % epoch start in ms
    S2.posttrig = 2000; % epoch end in ms   
    
    S2.trialdef(1).conditionlabel = 'Neutral face';
    S2.trialdef(1).eventtype = 'STI101_down';
    S2.trialdef(1).eventvalue = 1;
    S2.trialdef(2).conditionlabel = 'Happy face';
    S2.trialdef(2).eventtype = 'STI101_down';
    S2.trialdef(2).eventvalue = 2;
    S2.trialdef(3).conditionlabel = 'Fearful face';
    S2.trialdef(3).eventtype = 'STI101_down';
    S2.trialdef(3).eventvalue = 3;
    S2.trialdef(4).conditionlabel = 'Motorbike';
    S2.trialdef(4).eventtype = 'STI101_down';
    S2.trialdef(4).eventvalue = 4;
    
    S2.reviewtrials = 0;
    S2.save = 0;
    S2.epochinfo.padding = 0;
    S2.event = D_continuous.events;
    S2.fsample = D_continuous.fsample;
    S2.timeonset = D_continuous.timeonset;
    
    [epochinfo.trl, epochinfo.conditionlabels, S3] = spm_eeg_definetrial(S2);    
    
    %%%%
    % adjust timings to account for delay between trigger and visual display 
    if(0)
        timing_delay=0.030; % secs
        epochinfo_new=epochinfo;
        epochinfo_new.trl(:,1:2)=epochinfo.trl(:,1:2)+round(timing_delay/(1/D_continuous.fsample));
        epochinfo=epochinfo_new;
    end;
    
    %%%%
    % do epoching
    S3=[];
    S3.D = D_continuous;     
    S3.epochinfo=epochinfo;
    [D good_trial_start_times] = osl_epoch(S3);
            
end;

%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%
%% NOW DO VISUAL ARTEFACT REJECTION (REQUIRES MANUAL INPUT!)
% This runs a Fieldtrip interactive tool.
%
% - Pass over the first interactive figure as it is the EOG channel - so just
% press the "quit" button.
%
% This will bring up another interactive figure which will show the
% magnetometers. You can choose the metric to display - it is best to stick
% to the default, which is variance. This metric is then displayed for 
% the different trials (bottom left), the different channels (top
% right), and for the combination of the two (top left). You need to use
% this information to identify those trials and channels with high variance
% and remove them.  
%
% - Remove the worst channel (with highest variance) by drawing a box
% around it in the top right plot with the mouse. 
% - Now remove the trials with high variance by drawing a box
% around them in the bottom left plot.
% - Repeat this until you are happy that there are no more outliers.
%
% - Press "quit" and repeat the process for the gradiometers.

%set filenames used in following step
for subnum = 1:length(fif_files), % iterates over subjects
    spm_files{subnum}=[workingdir '/efA' spm_files_basenames{subnum}];
    %spm_files{subnum}=[workingdir '/efAd' spm_files_basenames{subnum}];
end

% RUN THE VISUAL ARTEFACT REJECTION:
for i=1:length(spm_file),
    S2=[];
    S2.D = spm_files{i};
    S2.time_range=[-0.2 0.4];
    D2=osl_rejectvisual(S2);
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% OPTIONAL EXTRAS BELOW:

%%%%%%%%%%%%%%%%%%%
%% Run this code to produce a plot of the number of rejected trials for each
% subject in the analysis
if(0)
    for i=1:length(spm_files),
        S2=[];
        spmfilename=spm_files{i};
        S2.D = spmfilename;   
        D=spm_eeg_load(S2.D);

        badchans(i)=sum(D.badchannels>0);
        rejected(i)=sum(D.reject);
    end;
    figure;plot(rejected,'*');xlabel('subject no.');ylabel('Num trials rejected'); % plot the number of rejected trials
    figure;plot(badchans,'*');xlabel('subject no.');ylabel('Num channels rejected'); % plot the number of rejected channels
end;
%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%
%% Run this bit of code if you want to unreject all trials and channels
if(0)
    for i= 1:length(spm_files),
        S2=[];
        spmfilename=spm_files{i};
        S2.D = spmfilename;        
        D=spm_eeg_load(S2.D);
 
        rejected=sum(D.reject)
        tmp=1:length(D.conditions);
        trlsel = zeros(1, length(tmp));
        D = reject(D, tmp, trlsel); 
        rejected=sum(D.reject)
              
        chans=D.meegchannels;
        D = badchannels(D,chans,zeros(size(chans)));
        save(D);
    end;
end;
%%%%%%%%%%%%%%%%%%%
